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Abstract: 

A grand canonical system of non-interacting fermions on a square lattice is considered at 
zero temperature. Three different phases exist: an empty lattice, a completely filled lattice 
and a liquid phase which interpolates between the other two phases. The Fermi statistics 
can be changed into a Bose statistics by coupling a statistical gauge field to the fermions. 
Using a quantum Monte Carlo method we investigate the effect of the gauge field on the 
critical properties of the lattice fermions. It turns out that there is no significant change 
of the phase diagram or the density of particles due to the gauge field even at the critical 
points. This result supports a recent conjecture by Huang and Wu that certain properties 
of a three-dimensional flux line system (which is equivalent to two-dimensional hard-core 
bosons) can be explained with non-interacting fermion models. 
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1. Introduction 

Non-interacting fermions play a fundamental role in the theory of exactly soluble models in 
one-dimensional quantum systems (e.g. hard-core bosons [1]) and two-dimensional classical 
systems (e.g. Ising spins [2], dimers [3] and directed interacting random walks [4]). It was 
suggested more recently that non-interacting fermions can also be used to model higher 
dimensional systems. Examples are magnets [5] and flux line systems of the mixed phase 
in high temperature superconductors with (unphysical) negative fugacity [6] in d = 3. On 
the other hand, a three-dimensional flux line system with (physical) positive fugacity is 
equivalent to a two-dimensional system of hard-core bosons [7]. Phase transitions (e.g. 
the Meissner-Abrikosov transition in a superconductor) exist for positive and for negative 
fugacities. Their properties are the same in two-dimensional flux line systems. This reflects 
the fact that one-dimensional hard-core bosons are equivalent to non-interacting fermions. 
Huang and Wu conjectured that the critical properties of thermodynamics are also the 
same [6] in three dimensions. For instance, the critical exponent for the density of flux 
lines n at the transition from the Abrikosov phase (n > 0) to the Meissner phase (n = 0) 
are the same. It is clear that certain properties of the hard-core bosons (e.g. superfluidity 
[8]) are not possible for non-interacting fermions. Therefore, the equivalence of flux line 
systems with negative and positive fugacity cannot be complete in three dimensions. 

The connection between fermions and bosons can be formally understood if we couple 
a statistical gauge field to the non-interacting fermions. The gauge field can change the 
Fermi statistics to Bose statistics. A well-known case is the Chern-Simons gauge field which 
transforms in two dimensions fermions into effective particles (anyons) [9] . The statistics of 
the anyons depends on the coupling constant. Another example for changing the statistics 
of fermions into bosons by means of a statistical field coupled to non-interacting fermions 
was discussed in Refs. 10-12. In that case two fermions (e.g., spin-dependent fermions with 
spin up and down) are bound in pairs by the gauge field. This construction is not restricted 
to two dimensions but can create hard-core bosons from fermions in any dimension. 

The effect of the statistical field was studied by means of large N methods in two extreme 
regimes: in the dense limit [10,12] and in the dilute limit [11]. Unfortunately, the prop- 
erties in the intermediate regime are not known. In particular, it would be interesting 
to understand whether there is a cross-over behavior from the dilute system to the dense 
regime or a transition between two different phases. For this purpose we performed a 
Monte Carlo simulation for the fermions in the statistical field. The method and some 
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results (density of particles and the phase diagram) are presented in the following. 

The article is organized as follows: The model of free lattice fermions and the coupling 
of a statistical gauge field are discussed in Sect. 2. Then we describe the quantum Monte 
Carlo method in Sect. 3. The numerical results for the density of hard-core bosons in 
d = 2 (or flux lines in d = 3) and the phase diagram are presented in Sect. 4. The results 
of the Monte Carlo simulation are compared in Sects. 5, 6 with the results of the free 
fermions of Sect. 2. 

2. The Model 

We define the dynamics of fermions on a square lattice with M sites by an inverse lattice 
propagator w which describes static fermions (resting particles) (t, r) — > (t + A, r) and 
hopping fermions (£, r) — > (t + A, r'). r and r' are the coordinates on the square lattice 
and t the (discrete) time. Subsequently we will also use space-time notation x = (t,r). 
The matrix elements of w are 

{1 for r' = r 
AJ for r', r nearest neighbors • (1) 
otherwise 

The partition sum, which is a generating function for physical quantities like the particle 
density, is the fermion determinant Z = det(w + (). ( is the fugacity for empty sites. Since 
we consider temperature T = 0, the free energy of the fermions reads as an integral over 
the Matsubara frequency to and the two-dimensional Brillouin zone [— 7r, n) 2 

/TT nTV nTV 
/ / log{e luJ [1 + t(cosfci + cos/c 2 )] + Qdk x dk 2 duj (2) 
-7T J —TT J —TT 

with i = 2AJ. Critical lines exist for wave vectors k\ = hi = (corresponding to a 
transition to an empty lattice) and k\ = k 2 = n (corresponding to a transition to a 
completely filled lattice) with ( = 1 ± 2t > 0. The density of fermions on a square lattice 
reads 

OF 

n = 1 -fac- (3) 

According to the definition of £ as the fugacity of empty sites the density is n = 1 for 
C = 0. Equations (2) and (3) lead, after some straightforward calculations, to 
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,"(C-i)/t 

n = l-— / 6(2 - |x|)K(l -x 2 /4)dx. (4) 

J-(C+i)/t 



6(2 — |x|) is the Heaviside step function and K(m) the elliptic integral [13]. The density n 
vanishes below t = (C - l)/2 (for ( > 1) and is n = 1 below t = (1 - 0/2 (for C < 1). This 
is consistent with the critical lines discussed above. Above the critical lines the density 
varies continuously between zero and one. This regime corresponds to a liquid state of 
fermions. The density of the fermions is shown as a function of 7 = y 7 ^ in Fig. (4). 

Usually a gauge field can be coupled to the lattice fermions by the transformation 

Wx,x> -»■ Ux,x' = e l * x ' x ' Wx,x' , (5) 

where the phase is anti-symmetric ($> x ,x' = —$>x',x) such that the field u XjX i is Hermitean. 
The fermions are transformed into hard-core bosons by choosing a real field u XjX > which 
is statistically independent for each pair of space-time neighboring sites x and x' [10-12]. 
Although the new field u is not Hermitean we will call it in this article a gauge field because 
it couples to the links of the lattice fermions. 

The partition sum of the hard-core bosons can be expressed as 

Z = (det 2 (u + 1 )) u . (6) 

The fermion determinant is squared because the gauge field is coupled to two independent 
fermions (e.g. spin 1/2 fermions). The fugacity is now 7 = ^/C- The statistical field has 
zero mean and the variance is 

( U t,r;t+A,r')u = W t ,r;t+A,r' ■ (7) 

The partition sum Z in Eq. (6) is the partition sum of flux lines or Bose world lines 
(BWLs). This becomes immediately clear if we expand the determinants and average with 
respect to the distribution of u: 

(det(-U + 7) 2 ) u = ^ (-^^(YllUx^iix) +7^x, 7 r 1 (x)]K,7r 2 (x) +l$ X ,n 2 (x)])u 

7Tl,7T2 X 

= (-ir 1 "^ 2 nt^x.Tr^x) + 7 2 ^x,7r 1 (x)]5 7ri ( x ),7r 2 (x) = y^TT[ W a,7ri(:c) + C^x,7n (x)] , (8) 

•7ri,7T2 X IT X 

where TTj(x) are permutations of space-time sites {x}. The last sum is indeed the partition 
sum of flux lines (or hard-core bosons) because it takes into account all possible configura- 
tion of paths from t = to t = (3 which are not intersecting. A flux line (or BWL) element 
has the weight w XjX '. (For details see Refs. 10-12.) 
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The density of BWLs at the inverse temperature (3 reads 

The factor in the partition sum Z can be canceled by rescaling Z — > Z = Z . 

Then £ _1 is the fugacity of bosons in Z , and the density reads with the effective chemical 
potential n = — log ( 

n = A^^^- (10) 
Eqs. (6), (7) and (10) are the starting points for our Monte Carlo simulation. 

3. The QMC - Method 

The form (6) of the partition sum allows us to perform a Monte Carlo - simulation of our 
system. This can be seen if one chooses as a (discrete) distribution P{u XjX i) a symmetric 
5 - function 

P(Ux,x') = l/2{5(Wx,x' + y/Wx,x') + &{ U x,x> - y/Wx,x')} (H) 

that satisfies the required properties (7). Any observable can be calculated by taking the 
average 

, T.\ u , l -4(w a; ,x',7)det 2 (w + 7) 
(A) = (A det 2 (u + 7 )) u = 7727-—^ • ( 12 ) 

For instance, to obtain the density of bosons 71(7, T) we have to set A = 1 — ^Tr[(u + 
r ))~ 1 ]/2J\f 13 according to Eq. (9). Of course, one cannot calculate Z as a function of u x x > 
and 7 for a reasonable lattice size. But since there are configurations of the u XjX i with 
very different weight with respect to their contribution to the partition sum, one should 
search for a way to sum only the main configurations, i.e., one should sum the weights of 
the configurations with respect to their probability to be realized. 

This is done by the Monte Carlo - algorithm. After choosing some initial configuration 
{u x ,x'} w e generate a Markov process by flipping elements u XjX > with respect to their 
former value with some probability W({u XyX >} — > {u x , x '}). This probability has to fulfill 
the property of 'detailed balance': 

W({u XiX/ } -> {«x/}>(K,x'}) = W({u XiX/ } -> (13) 

where p({u XiX /}) denotes the equilibrium probability of the corresponding configuration 
{u XjX '}- We used the minimum of some random number of a rectangular distribution be- 
tween and 1 and the 'acceptance ratio' R (the squared ratio of the determinants of the new 
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and the old configurations of u XjX i) as the transition probability (Metropolis algorithm). 
The acceptance ratio R can be calculated via Green's function, the inverse of the matrix 
(u + ~i) x ,x'- If flipping of an element u x , x > is accepted we can update Green's function 
with low computational effort. We are repeating this procedure ('sweeping' through the 
space-time lattice) until we are sure to be close to the equilibrium configuration. Then we 
can begin to measure the desired quantity and repeat this as long as necessary to produce 
stable results with low statistical error. 

This QMC - algorithm is related to the usual fermionic QMC - method established by 
Blankenbecler et al. [14,15]. Since the acceptance ratio is a square, it is always positive 
and the fermionic sign problem does not appear, equivalent, e.g., to the fermionic Hubbard 
model at half filling. The main problem of this QMC - method lies in the fact that the 
flipping variable fixed to the bonds of the lattice and not to the sites. Obviously, 

the number of flipping variables is (2d + l)Af, since a particle can either hop (2d) or stay at 
its spatial point (1), but it has to go on in the time direction. This restricts the simulations 
to small lattice sizes (max. M = 100 sites). 

The procedure has been checked by comparing the kinetic energy at half density with the 
exact solution (d = 1) [16], exact diagonalization and QMC - data [17,18] of the quantum 
spin- 1/2 XY- model to which the hard-core bosons can be mapped [19]. We find good 
agreement which establishes the validity of our algorithm. 

4. Numerical results 

i) Phase transitions at J = 1 

The numerical work was mostly done on the Siemens VP-600EX of the RZ Karlsruhe. A 
typical job for the 8x8 lattice at low temperatures needs about one hour CPU time. The 
results shown are for a simple square lattice in dimension d = 2. We also simulated bosons 
on a triangular lattice. After adjusting the hopping strength we find that the system 
behaves as in the case of a square lattice. This was expected since there is no frustration 
in both cases. 

First we consider the dependence of the density on the fugacity 7 which is the natural 
parameter for the simulation and from which we can deduce the dependence of the density 
on the bosonic chemical potential \i. We expect from the results in [10,12] that a superfluid- 
insulator transition takes place at the critical value 7 C = 1 at T = in an infinite system. 
Due to finite lattice sizes and finite temperature T in our QMC-simulation the sharp 
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transition is rounded. To estimate whether a sharp phase transition really exists we must 
extrapolate our data to T = and M — > oo. This is shown in Figs. 1 and 2 . The curves 
for 7 < 7 C saturate at a finite density below a 7-dependent temperature (Fig. 1). In 
contrast, the graph for 7 C does not saturate but decreases monotonically with dominant 
linear behavior which can be extrapolated to a finite density at T = 0. This behavior is 
characteristic for 7 C for a given hopping strength J and can be used to determine the phase 
boundaries (see below). The graphs for 7 > 7 C (only 7 = 1.005 is shown) decrease faster 
and the density becomes exponentially small as T — > even for finite lattice sizes. 
To analyze finite size effects we look at the density for different lattice sizes M (Fig. 2). 
As expected, the density generally decreases upon increasing the lattice size. The density 
for 7 = .995 < 7 C is clearly finite in the thermodynamic limit {M — > 00), even at T = 0. 
However, the remaining density at T = for the critical value 7 C = 1 scales to zero as 1/Af 
(see inset of Fig. 2) within numerical accuracy. 

We therefore conclude that the numerical data indicate a second order phase transition at 
T = in the thermodynamic limit for J = 1 at 7 C = 1, in accord with analytical results 
[10,12]. 

ii) Phase diagram 

To obtain the phase diagram in the 7 - J - plane we have to vary J and look for the value 
of 7 at which the density is linearly decreasing at low temperatures but extrapolates to a 
finite remaining density at T = 0. Since we do not know that value a priori we have to do a 
considerable amount of jobs to estimate accurate data. Due to limited CPU-time available 
we restrict ourselves to the 4x4 lattice to calculate the shape of the phase boundaries and 
then check some points in the 8x8 lattice. We find no significant changes of the boundaries 
in the bigger system. 

First we observe a monotonic shrinking of the superfluid region as the hopping J is de- 
creased as can be seen from Fig. 3. At hopping J = 0, we know that there is a first order 
phase transition between the commensurate phases with n = and n = 1. For very small 
hopping (J <.02) the superfluid region becomes too small to distinguish the numerical 
data (with statistical error less than 2%) from the exact data at zero hopping even for 
very low temperatures. A possible phase diagram might show an extension of this first or- 
der point to values of J > 0, forming a first order transition line separating the two phases 
commensurate with the lattice. The numerical data do not rule out such an extension, 
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but if it exists it must be very small (J m ax < -02). We assume in the reasoning below 
that there is no extension, but the expressions can be easily modified if this is not valid. 
Therefore, we obtain the phase diagram shown as Fig. 3. The Monte Carlo data (circles) 
fit well with the phase boundaries for a system with noninteracting fermions. Especially 
the transitions from the empty lattice to a finite density are in very good agreement with 
each other, even for much larger values of J than shown. The transitions from finite den- 
sity to an insulating state with density n = 1 are in good agreement for moderate J, but 
deviations become visible at about J = 1. Not only the phase boundaries, but also the 
density profiles n vs. 7 of hard-core bosons and free fermions are in good agreement for 
small hopping, see Fig. (4). For larger hopping J there are deviations, e.g. the curvature 
of the boson and fermion data is different for intermediate densities. Still, the qualitative 
agreement is quite good. 

We found a second order transition at finite J but a density jump of the size of unity for 
J = 0. This suggests a scaling ansatz for the density 

noc( 7 ~ 7c(J V (14) 

rJ 

with rj = 1 to fit the second order transition just to the right of the first order point at 
J c = 0. In the vicinity of J c the behavior (14) should be dominant. Indeed, we observe 
a strong increase of the (negative) slope of the density profile for decreasing J (Fig. 5). 
If we plot the slope from the extrapolation of our data to T = 0, we find indeed a linear 
dependence of the slope on the inverse hopping 1/J (see inset of Fig. 5) for moderate 
values of J. This is in agreement with the ansatz (14). 

5. Conclusions 

The numerical results give a complete picture of the phase transitions in a pure grand 
canonical system of three-dimensional flux lines (or two-dimensional hard-core bosons) on 
a lattice. We found no sign for a phase which is related to a special gauge field configuration 
(e.g., a "flux phase" found for the t — J model [20]). The data indicate second order phase 
transitions at nonzero hopping J (with the caveat about the first order line for very small 
hopping) from an empty system (density n = 0) to a superfluid phase with noninteger n 
and finally to an insulating phase with density n = 1 as we increase the chemical potential 
(decrease ( or 7). This behavior was also found for the free fermion system given in 
Eq. (4). The 7 — J - phase diagram of the free fermions agrees with that of the hard- 
core bosons within numerical errors (Fig. (3)). No sign for a singular behavior of the 
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boson system between the dense and the dilute regime was observed. A reasonably good 
agreement of fermion and boson systems is also found for the density profile (c.f. Fig. 
4). Therefore, we conclude that the thermodynamic properties (i.e., quantities which can 
be derived from (2) and (6), respectively) are not affected by the gauge field fluctuations 
u. This supports the conjecture of Huang and Wu that the critical properties of higher 
dimensional systems (here: three-dimensional flux lines and two-dimensional hard-core 
bosons) can be described by a free fermion statistics. This relation offers a new approach 
for the statistics of flux-lines (or hard-core bosons) in random potential where we expect 
new physical states like the vortex glass [21] and the Bose glass [22,23]. For instance, 
the freezing of the fermionic line dynamics due to disorder [24] might be related to the 
vortex glass transition of physical flux lines. The analogy of free fermions and flux lines 
for d > 2 would also allow the extension of the two-dimensional results, where the creation 
of overhangs and finite loops of flux lines by replica symmetry breaking was found [25]. 
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Figure Captions 

Fig. 1: Density of bosons versus temperature for various fugacities 7onaA/" = 8x8 = 64 
lattice (hopping J = 1). The density saturates for 7 < 7 C = 1 and becomes exponentially 
small for 7 > 7 C as temperature T — > 0. For 7 = 7 C the density drops linearly and 
approaches a finite, lattice size dependent value. This behavior is characteristic and allows 
us to find the critical 7 C for arbitrary hopping strength J. T Q is a reference temperature. 

Fig. 2: Density of bosons versus temperature for various lattice sizes Af for the fugacity 
7 = 0.995 and 7=1, respectively (hopping strength J = 1). The density is a decreasing 
function of lattice size for all temperatures. Extrapolating to T = gives density values 
which stay finite as Af is growing for 7 < 7 C . For 7 C = 1 the density at T = behaves 
like 1/Af and extrapolates to zero as Af — > 00 (see inset). This supports the picture of a 
second order phase transition at T = and 7 C = 1. T a is a reference temperature. 

Fig. 3: Phase diagram of hard-core bosons and free fermions on a square lattice in the 
7 — J plane. The circles represent the QMC data for hard-core bosons, the dashed lines 
are the phase boundaries for noninteracting, spinless fermions. The phases are insulating 
(Mott insulator) above the upper (density n = 0) and below the lower phase boundary 
(n = 1). In between, the particles are in a (super-) fluid phase. The agreement of the 
phase boundaries of the two systems is good for small hoppings J and also for large J for 
the upper boundaries. The lower boundaries deviate for larger J. 

Fig. 4: Density profile for various hoppings J for hard-core bosons and free fermions on a 
square lattice. One clearly observes the (insulating) phases with density n = and n = 1, 
respectively. There is qualitative agreement of bosons (symbols) and fermions (solid lines) 
for small hoppings but obvious deviations for larger hopping J. The temperatures are in 
units of the reference temperature T Q . 

Fig. 5: Density versus fugacity 7 for three intermediate values of hopping J. The slope 
—dn/d'y increases upon decreasing J. The inset shows that the slope is roughly propor- 
tional to the inverse hopping 1/J. The value of the slopes in the inset are determined 
from the extrapolation of the data to T = 0. The dashed lines are a guide to the eye. The 
temperature is T = 1/6T G . 
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